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Abstract In this paper, we present a multigrid preconditioner for solving 

^~~^ the linear system arising from the piecewise linear nonconforming Crouzeix- 

Raviart discretization of second order elliptic problems with jump coefficients. 

I— —I The preconditioner uses the standard conforming subspaces as coarse spaces. 

-^ Numerical tests show both robustness with respect to the jump in the coeffi- 

^^ cient and near-optimality with respect to the number of degrees of freedom. 

ri 1 Introduction 

, The purpose of this paper is to present a multigrid preconditioner for solving 

^ the linear system arising from the P^ nonconforming Crouzeix-Raviart (CR) 

f^ discretization of second order elliptic problems with jump coefficients. The 

^<0 multigrid preconditioner we consider here uses pointwise relaxation (point 

'^ Gauss-Seidel/Jacobi iterative methods) as a smoother, followed by a subspace 

. (coarse grid) correction which uses the standard multilevel structure for the 

C^^ nested Pi conforming finite element spaces. The subspace correction step is 

^— ^ motivated by the observation that the standard P^ conforming space is a 

I subspace of the CR finite element space. 

• • One of the main benefits of this algorithm is that it is very easy to ini- 

. J^ plement in practice. The procedure is the same as the standard multigrid 

S^ algorithm on conforming spaces, and the only difference is the prolongation 

\—t and restriction matrices on the finest level. Since the spaces are nested, the 
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prolongation matrix is simply the matrix representation of the natural inclu- 
sion operator from the conforming space to the CR space. 

The idea of using conforming subspaces to construct preconditioners for 
CR discretization has been used in Xu [1989, 1996] in the context of smooth 
coefficients. For the case of jumps in the coefficients, domain decomposition 
preconditioners have been studied in Sarkis [1994b, a] and the BPX precondi- 
tioner has been considered in Ayuso de Dios et al. [2010] in connection with 
preconditioners for discontinuous Galcrkin methods. 

In the context of jump coefficients, the analysis of multigrid precondition- 
ers for conforming discretizations is given in Xu and Zhu [2008]. For CR 
discretizations, the analysis is more involved due to the nonconformity of the 
space, and special technical tools developed in Ayuso de Dios et al. [2010] are 
necessary. Due to space restrictions, we only state the main result (Theorem 2 
in Section 3), and provide numerical results that support it. Detailed analyses 
and further discussion of the algorithm will be presented in a forthcoming 
paper. 

The paper is organized as follows. In Section 2, we give basic notation and 
the finite element discretizations. In Section 3, we present the multigrid algo- 
rithm and discuss its implementation and convergence. Finally, in Section 4 
we verify numerically the theoretical results by presenting several numerical 
tests for two and three dimensional model problems. 



2 Preliminaries 

Let n C R"^ {d — 2,3) be an open polygonal domain. Given / G L^(J?), we 
consider the following model problem: Find u E Hl{fi) such that 

a{u,v):^{nVu,Vv) = {f,v) ^v e H^i^) , (1) 

where the diffusion coefficient k € L°°{fi) is assumed to be piecewise con- 
stant, namely, K,{x)\a^ = Km is a constant for each (open) polygonal subdo- 
main i?™ satisfying Um^i^m — ^ and ilm H i7„ == for m 7^ n. 

We assume that there is an initial (quasi-uniform) triangulation To, with 
mesh size Hq, such that for all T G To kt '■— '*(2;)|t is constant. Let T) '■— Th 
(j = 1, • • • , J) be a family of uniform refinement of To with mesh size hj. 
Without loss of generality, we assume that the mesh size hj ~ 2^^ ho (j = 
0, • • • , J) and denote h = hj. 

On each level j — 0, • • • , J, we define Vj as the standard P^ conform- 
ing finite element space defined on T). Then the standard conforming finite 
element discretization of (1) reads: 

Find Uj G Vj such that a{uj,Vj) — {f,Vj), Vwj G Vj. (2) 

For each j = 0, • • • , J, wc define the induced operator for (2) as 
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{AjVj,wj) = a{vj,Wj), yv.j,Wj e Vj. 

We denote £h the set of all edges (in 2D) or faces (in 3D) of Th- Let V^^ 
be the piecewise linear nonconforming Crouzeix-Raviart finite element space 
defined by: 



vr- 



■IveL^in) : v\^ e ¥\T)yTeTh and [{vjeds = OW e e £h\ , 



where ¥^{T) denotes the space of linear polynomials on T and |w|e denotes 
the jump across the edge/face e € £h with \v\e = v when e C dil. In the 
sequel, let us denote Vj+i := V^^ for simplicity. We remark that all these 
finite element spaces are nested, that is, 

VoC--- cVj c Vj+1. 

The P^-nonconforming finite element approximation to (1) reads: 

Find u e V^^ : ah{u,w) -^ J2 [ '*^^" ' '^"' = {.f,w),yw € V,p^. (3) 

TeTj '^ 

The bilinear form a/j(-, •) induced a natural energy norm: \v\h^K, '■= \/cih{v, v) 
for any v G V^'^. In operator form, we are going to solve the linear system 

Au = /, (4) 

where A is the operator induced by (3), namely 

{Av,w) ^ ah{v,w), yv,w eV,^'^. 



3 A Multigrid Preconditioner 

The action of the standard multigrid V-cycle preconditioner B := -Bj+i : 
Vj+i I— > Vj+i on a given g G Vj+i is recursively defined by the following 
algorithm (cf. Bramble [1993]): 

Algorithm 1 (y-cycle) Let g,j+i — g, and B^ ~ A^^ . For j = 1, • • • , J+1, 
we define recursively Bjgj for any gj G Vj by the following three steps: 

1. Pre-smoothing : wi = Rjgj] 

2. Subspace correction: 1/72=^1+ Bj-iQj^i{gj — AjWi); 

3. Post-smoothing: Bjgj :~ W2 + R*j{9j ^ AjW2). 

In this algorithm, Rj corresponds to a Gauss-Seidel or a Jacobi iterative 
method known as a smoother; and Qj is the standard L^ projection on Vj: 
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(QjV, Wj) = (v, Wj), VWj eVj, (j = 0, • • • , J). 

The implementation of Algorithm 1 is almost identical to the implementa- 
tion of the standard multigrid y-cyclc (cf. Briggs et al. [2000]). Between the 
conforming spaces, we use the standard prolongation and restriction matrices 
(for conforming finite elements) . The corresponding matrices between Vj and 
Vj+i, arc however different. The prolongation matrix on Vj can be viewed 
as the matrix representation of the natural inclusion Xj : Vj ^i- Vj^i, which 
is defined by 

(Xju)(x) = ^ V{me)ll}e{x), 
eeSh 

where ij^e is the CR basis on the edge/face e G £h and nie is the barycenter 
of e. Therefore, the prolongation matrix has the same sparsity pattern as 
the edge- to- vertex (in 2D), or face- to- vertex (in 3D) connectivity, and each 
nonzero entry in this matrix equals the constant 1/d where d is the space 
dimension. The restriction matrix is simply the transpose of the prolongation 
matrix. 

The efficiency and robustness of this preconditioncr can be analyzed in 
terms of the effective condition number (cf. Xu and Zhu [2008]) defined as 
follows: 

Definition 1. Let F be a real N dimensional Hilbert space, and S : V —> V 
be a symmetric positive definition operator with eigenvalues < Ai < • • • < 
Xn ■ The TO-th effective condition number of S is defined by 

Note that the standard condition number IC{BA) of the preconditioned sys- 
tem BA will be large due to the large jump in the coefficient k. However, 
there might be only a small (fixed) number of small eigenvalues of BA, which 
cause the large condition number; and the other eigenvalues are bounded 
nearly uniformly. In particular, we have the following main result: 

Theorem 2. Let B be the multigrid V-cycle preconditioncr defined in Algo- 
rithm 1. Then there exists a fixed integer toq < M, depending only on the 
distribution of the coefficient k, such that 

/C™„(BA)<C2|log/i|2 = cV2, 

where the constant C > is independent of the coefficients and mesh size. 

The analysis is based on the subspace correction framework Xu [1992], but 
some technical tools developed in Ayuso de Dios et al. [2010] are needed to 
deal with nonconformity of the finite element spaces. Due to space restriction, 
a detailed analysis will be reported somewhere else. 

Thanks to Theorem 2 and a standard PCG convergence result (cf. [Ax- 
elsson, 1994, Section 13.2]), the PCG algorithm with the multigrid F-cycle 
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preconditioncr defined in Algorithm 1 has the following convergence estimate: 

\u - «.U,« < 2{UBA) - l)™" [-^j^A I" - «ok.« , 

where Mq is the initial guess, and Ui is the solution of i-th PCG iteration. 
Although the condition number IC{BA) might be large, the convergence rate 
of the PCG algorithm is asymptotically dominated by §jxy, which is deter- 
mined by the effective condition number ICmoiBA). Moreover, this bound of 
asymptotic convergence rate convergence is independent of the coeffieient k, 
but depends on the mesh size logarithmically. 



4 Numerical Results 

In this section, we present several numerical tests in 2D and 3D which verify 
the result in Theorem 2 on the performance of the multigrid V^-cycle precon- 
ditioncr described in the previous sections. The numerical tests show that 
the effective condition numbers of the preconditioned linear systems (with 
F-eycle preconditioner) are nearly uniformly bounded. 



4-1 A 2D Example 

As a first model problem, we consider equation (1) in the square f2 = (—1, 1)^ 
with coefficient such that, k{x) = 1 for x e i7i = (—0.5,0)^ U (0,0.5)^, and 
k{x) = e for X in the remaining subdomain, x G J? \ J7i (see Figure 1). By 
decreasing the value of e we increase the contrast in the PDE coefficients. 

Our initial triangulation on level has mesh size ho ~ 2^^ and resolves the 
interfaces where the coefficients have discontinuities. Then on each level, we 
uniformly refine the mesh by subdividing each element into four congruent 
children. In this example, we use 1 forward/backward Gauss-Seidel iteration 
as pre/post smoother in the multigrid preconditioner, and the stopping cri- 
teria of the PCG algorithm is ||?'fc||/||^o|| < 10^^ where r^ is the the residual 
at fc-th iteration. 

Figure 2 shows the eigenvalue distribution of the multigrid V^-cycle pre- 
conditioned system BA when h = 2~^ (level =4) and e — 10^^. As we can 
see from this figure, there is only one small eigenvalue that deteriorates with 
respect to the jump in the coefficient and the mesh size. 

Table 4.1 shows the estimated condition number /C and the effective con- 
dition number /Ci of BA. It can be observed that the condition number /C 
increases rapidly with respect to the increase of the jump in the coefficients 
and the number of degrees of freedom. On the other hand, the number of PCG 
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Fig. 1 2D Computational Domain 



Fig. 2 Eigenvalue Distribution of BA 



iterations increases only a small amount, and the corresponding effective con- 
dition number is nearly uniformly bounded, as predicted by Theorem 2. 
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levels 





1 


2 


3 


4 


1 




1.65 (8) 
1.44 


1.83 (10) 
1.78 


1.9 (10) 
1.77 


1.9 (10) 
1.78 


1.89 (10) 
1.76 


10-1 


K 
ICi 


3.78 (10) 
1.89 


3.69 (11) 
1.87 


3.76 (12) 
1.93 


3.79 (12) 
1.92 


3.88 (12) 
1.95 


10-2 


K 


23.4 (12) 
2.15 


23.6 (13) 
1.96 


24.6 (13) 
1.99 


25.1 (14) 
1.97 


26 (15) 
2.24 


10-3 


K. 
/Ci 


218 (13) 
2.19 


223 (14) 
1.98 


232 (15) 
2 


238 (16) 
1.98 


246 (16) 
2.29 


10-4 


K 


2.17e-|-03 (14) 
2.2 


2.21c-(-03 (15) 
1.98 


2.31c-(-03 (16) 
2 


2.37e-|-03 (18) 
1.98 


2.45e+03 (18) 
2.3 


10-5 




2.17C-I-04 (15) 
2.2 


2.21e-|-04 (16) 
1.98 


2.31e-|-04 (17) 
2 


2.37e-|-04 (19) 
1.98 


2.76e+04 (19) 
2.64 



Table 1 Estimated condition number K (number of PCG iterations) and the effective 
condition number fCi 



4.2 A 3D Example 



In this second example, we consider the model problem (1) in the open unit 
cube in 3D with a similar setting for the coefficient. We set k{x) = 1 for 
X £ f2i — (0.25, 0.5)'^ or a; G i72 = (0.5, 0.75)^, and k{x) — e for the remaining 
subdomain (that is, for x € f2\(f2iUf22))- The domain i? and the subdomains 
just described are shown in Figure 3. The coarsest partition has mesh size 
ho = 2-^, and it is set in a way so that it resolves the interfaces where the 
coefficient has jumps. 
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To test the effects of tlie smoother, in this example we used 5 for- 
ward/backward Gauss-Seidel as smoother in the multigrid preconditioner. 
In order to test more severe jumps in the coefficients, we set the stopping 
criteria ||?'fe||/||''o|| < 10~^^ for the PCG algorithm in this experiment. 




^1 




Fig. 3 3D Computational Domain 



Fig. 4 Eigenvalue Distribution of BA 



Figure 4 shows the eigenvalue distribution of the multigrid 1^-cycle precon- 
ditioned system BA when h — 2^^ (level=3) and e = 10~^. As before, this 
figure shows that there is only one small eigenvalue that even deteriorates 
with respect to the jump in the coefhcients and the mesh size. 
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1.19 (8) 
1.16 


1.34 (11) 
1.26 


1.37 (11) 
1.31 


1.36 (11) 
1.29 


10-1 
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2.3 (10) 
1.60 


1.94(13) 
1.56 


1.75 (13) 
1.45 


1.67 (14) 
1.43 


10-^ 


K. 


86.01 (11) 
2.4 


63.07 (16) 
2.12 


52.67 (17) 
1.89 


48.19(17) 
1.78 


10-5 
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8.39+03 (13) 

2.44 


6.15e+03 (18) 
2.14 


5.13C+03 (19) 
1.91 


4.70e+03(19) 
1.80 


10-7 


K 
^1 


8.39+05 (14) 
2.45 


6.15e+05 (21) 
2.14 


5.13e+05 (23) 
1.91 


4.70e+05(21) 
1.80 



Table 2 Estimated condition number K (number of PCG iterations) and effective condi- 
tion number K\. 



Table 2 shows the estimated condition number /C (with the number of PCG 
iterations), and the effective condition number K-i. As is easily seen from the 
results in this table, the condition number /C increases when e decreases, i.e. 
the condition number grows when the jump in the coefficients becomes larger. 
On the other hand, the results in Table 2 show that the effective condition 
number K-i remains nearly uniformly bounded with respect to the mesh size 
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and it is robust with respect to the jump in the coefficient, thus confirming 
the result stated in Theorem 2: a PCG with multigrid T^-cycle preconditioner 
provides a robust, nearly optimal solver for the CR approximation to (3). 
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